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Abstract 



The exact A^-particle ground state wave function for a one-dimensional con- 
densate of hard core bosons in a harmonic trap is employed to obtain accurate 
numerical results for the one-particle density matrix, occupation number dis- 
tribution of the natural orbitals, and momentum distribution. Our results 
show that the occupation of the lowest orbital varies as 

^0.59 

in contrast to 

^0.5 £qj, ^ spatially uniform system, and for a true BEG. 
03.75.Fi,03.75.-b,05.30.Jp 



1 



Typeset using REVT^ 



I. INTRODUCTION 



Recent advances in atom de Broglie waveguide technology and its potential applica- 
bility to atom interferometry p and integrated atom optics create a need for accurate 
theoretical modelling of such systems in the low temperature, tight waveguide regime where 
transverse excitations are frozen out and the quantum dynamics becomes essentially one- 
dimensional (ID). It has been shown by Olshanii f^, and also recently by Petrov et al. 
0, that at sufficiently low temperatures, densities and large positive scattering length, a 
Bose-Einstein condensate (BEC) in a thin cigar-shaped trap has dynamics which approach 
those of a ID gas of hard core, or impenetrable, point bosons. This is a model for which the 
exact many-body energy eigensolutions were found in 1960 using an exact mapping from 
the Hilbert space of energy eigenstates of an ideal gas of fictitious spinless fermions to that 



of many-body eigenstates of hard core, and therefore strongly interacting, bosons pO| , p 
In this limit there are strong short-range pair correlations which are omitted in the Gross- 
Pitaevskii (GP) approximation. In the absence of a trap potential it is known [|T2| that the 
occupation of the lowest orbital is of order where N is the total number of atoms, in 
contrast to for the ideal Bose gas and GP approximation. Nevertheless, this system ex- 
hibits some BEC-like behavior such as Talbot recurrences following an optical lattice pulse 
||TB| and dark soliton-like behavior in response to a phase-imprinting pulse . 

The case of harmonically trapped, hard core bosons in ID is more relevant to recent 
atom waveguide experiments [|T3|. The spatial density profile of the single-particle density 
is expressible in closed form, and has recently been shown [|16| to be well approximated 
by a modified ID effective field theory, although we have recently shown in a numerically 
accurate time-dependent calculation that spatial interference of separated and recom- 
bined condensates is much weaker than that predicted by the corresponding time- dependent 
mean field theory [|1^]. Although the Fermi-Bose mapping theorem |]10|,|11[ implies that all 
physical properties expressible in terms of spatial configurational probabilities are the same 
for the actual bosonic system and the fictitious "spinless fermion" system used for the map- 
ping, the momentum distribution of the bosonic system, or more generally its occupation 
distribution over the relevant orbitals for a given geometry, is very different in the bosonic 



system. It is known p,12,18| that for a spatially uniform system of hard core bosons in 



ID, the momentum distribution is strongly peaked in the neighborhood of zero momentum, 
whereas that of the corresponding Fermi system is merely a filled Fermi sea. In the case 
of hard core bosons in a ID harmonic trap, it is an interesting and previously unanswered 
question whether the system undergoes true BEC or merely an attenuated one such as that 
in the uniform system. Ketterle and Van Druten have shown that true BEC occurs for 



a finite number of atoms in a ID harmonic oscillator (HO) for an ideal gas. We examine 
the question for hard core bosons in a ID HO by using the Fermi-Bose mapping theorem 
to generate the exact many-body ground state. The most fundamental definition of BEC 
and the condensate orbital is based on the large distance behavior of the one-particle re- 
duced density matrix pi{x,x'). If off-diagonal long-range order (ODLRO) is present and 
hence the largest eigenvalue of pi is macroscopic (proportional to A^) then the system is 
said to exhibit true BEC and the corresponding eigenfunction, the condensate orbital, plays 
the role of an order parameter [^,^. Although the precise definition of ODLRO requires 
a thermodynamic limit not strictly applicable to mesoscopic traps, the GP approximation 
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assumes from the start that ODLRO and macroscopic occupation of a single orbital are 
good approximations in a trap, so examination of this assumption in the Olshanii limit P] is 
important. In the remainder of this paper we shall determine the many-body ground state 
and its salient features, including the one-particle reduced density matrix and its eigenvalues 
(occupation number distribution function) and eigenfunctions (natural orbitals), as well as 
the momentum distribution function. 



II. EXACT GROUND STATE WAVE FUNCTION 



The Hamiltonian of N bosons in a ID harmonic trap is 



N 



1 2 2 



(1) 



We assume that the two-body interaction potential consists only of a hard core of ID diam- 
eter a. This is conveniently treated as a constraint on allowed wave functions ipi^ir ' ' ^^n) 
such that 

ip = if \xj-Xk\ <a , I < j < k < N, (2) 

rather than as an infinite interaction potential. It follows from the Fermi-Bose mapping 
theorem ||lO|JTT| , p^ that the exact N-boson ground state ipBo of the Hamiltonian (1) with the 
constraint (2) is 

ipBoixi, ■ ■ ■ ,xn) = \i^Fo{xi, ■ ■ ■,xn)\, (3) 

where ippo is the ground state of a fictitious system of N spinless fermions with the same 
Hamiltonian (|l]) and constraint. At low densities it is sufficient [^,^ to consider the case 
of impenetrable point particles, the zero-range limit a -h> of (2). Since wave functions of 
"spinless fermions" are antisymmetric under coordinate exchanges, their wave functions van- 
ish automatically whenever any xj = x^, the constraint has no effect, and the corresponding 
fermionic ground state is the ground state of the ideal gas of fermions, a Slater determinant 
of the lowest N single-particle eigenfunctions 0„ of the harmonic oscillator (HO) 

I {N-1,N) 

V'Fo(a;i,---,a;7v) = det M^j)- (4) 

\/Nl {nj)=(0,l) 



The HO orbitals are 



M^) = -^^7^-=e-«^/^/7„(g) (5) 



with HniQ) the Hermite polynomials and Q = x/xosc^osc = \J'h/muj being the ground 
state width of the harmonic trap for a single atom. By factoring the Gaussians out of the 
determinant and carrying out elementary row and column operations, one can cancel all 



terms in each except the one of highest degree, with the result |22 
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{N-1,N) {N-1,N) 

det HJxA = 2^(^-i)/2 det (x,)" 

(n,j)={0,l) ^ (n,i)={0,l) 



2N{N-l)/2 



(x^ Xj^ 

^<j<k<N 



(6) 



Substitution into (3) then yields a simple but exact analytical expression of Bijl-Jastrow 
pair product form for the A^-boson ground state: 



1) 



N 



■ N 

.1=1 



2/2 



n 

^<j<k<N 



X}^ X 



(7) 



with Qi = Xi/xosc and normalization constant 

1 



C 



N 



2N(N~l)/A 



Xn 



N-1 



,-1/2 



n=0 



(8) 



It is interesting to note the strong similarity between this exact ID A^-boson wave function 
and the famous Laughlin variational wave function of the 2D ground state for the quantized 
fractional Hall effect p3| , as well as the closely-related wave functions for bosons with weak 
repulsive delta-function interactions in a harmonic trap in 2D found recently by Smith and 
Wilkin [El. 



III. GROUND STATE PROPERTIES 

In this section we numerically evaluate the ground state properties of a ID condensate 
of N hard core bosons in a harmonic trap using the exact many-body wave function of the 
previous section. 



A. Single particle density and pair distribution function 

Both the single particle density and pair distribution function depend only on the ab- 
solute square of the many-body wave function, and since iV^soP = I'^foP they reduce to 
standard ideal Fermi gas expressions. The single particle density, normalized to A^, is 



p{x) = N lijjBoix, X2, ■ ■ ■ , XN)\^dx2 ■ ■ ■ dx 



N 



N-1 



E \Vn{x)\' (9) 



n=0 



We shall not exhibit it here since it has recently been calculated by Kolomeisky et al. [16 



see also our recent discussion of the time-dependent case [JT^. The pair distribution function, 
normalized to N[N — 1), is 



D{xi,X2) = N{N -I) I \iJBo{xi,---,XN)\ 



'^dx'x ■ ■ ■ dx 



N 



\(Pn{Xl)'^n'{x2) - (Pn{x2)'^n'{Xi)\^ (lO) 
0<n<ra'<Af-l 
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Physically, the pair distribution function is the joint probability density that if one atom is 
measured at Xi then a second measurement immediately following the first finds an atom at 
X2- Noting that terms with n = n', which vanish by antisymmetry, can be formally added 
to the summation (9), one can rewrite the pair distribution function in terms of the single 
particle density and a correlation function A: 



|2 



D{xi,X2) = p{xi)p{x2) - \ A{xi,X2) 

A{xi,X2) = f*niXl)fn{x2) (H) 



n=0 



Although the Hermite polynomials have disappeared from the expression (7) for the many- 
body wave function, they reappear upon integrating iV'soP over (A^ — 1) coordinates to 
get the single particle density p{x) and over (iV — 2) to get the pair distribution function 
D{xi,X2), and the expressions in terms of the HO orbitals (pn are the most convenient for 
evaluation. 

Figure shows a gray scale plot of the dimensionless pair distribution function xl^^ ■ 
D{Qi,Q2) versus the normalized coordinates Q12 = Xi^2/xosc for a) = 2, b) = 6, and 
c) N = 10. Some qualitative features of the pair distribution function are apparent: In 
the first place it follows either from the original expression (9) or from Eqs. (8) and (10) 
that D{xi,X2) vanishes at contact xi = X2, as it must because of impenetrability of the 
particles, and we see this to be true in Fig. |I[ Furthermore, the correlation term A{xi,X2) 
is a truncated closure sum and approaches the Dirac delta function 6{xi — X2) as N 00, as 
is to be expected since the healing length in a spatially uniform ID hard core Bose gas varies 



inversely with particle number As a result the width of the null around the diagonal 
Qi = Q2 decreases with increasing A^, and vanishes in the limit. Away from the diagonal 
along Q2 = —Qi the pair distribution function rises, exhibits modulations for N > 2, due to 
the oscillatory nature of the HO orbitals, before decreasing back to zero at large distances. 
For \xi — X2\ much larger than the healing length, D reduces to the uncorrelated density 
product p{xi)p{x2), so the spatial extent of the pair distribution function is that of the 
density and varies as N^^"^ |1^ . 



B. Reduced single-particle density matrix 

The reduced single-particle density matrix with normalization / pi{x,x)dx = N is given 



by 



Pl{x, x') = N J 1pBo{x, X2, ■ ■ ■ , Xn) 

xipBo{x', X2, ■ ■ ■ , Xj\f)dx2 ■ ■ ■ dx]\r 



r ^ 



i=2 

x[ n {Qk-Q,fYQ2---dQ^, (12) 



2<i<fc<A^ 

with 
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N 



Af-l 



n=0 



(13) 



Although the multi-dimensional integral (|T2D cannot be evaluated analytically, it can be eval- 
uated numerically by Monte Carlo integration for not too large values of N (the computing 
time scales as A^"^). Figure |^ shows a gray scale plot of the dimensionless reduced single- 
particle density matrix Xosc ■ Pi{.Q^Q') versus the normalized coordinates Q and Q' for a) 
AT = 2, b) A^ = 6, and c) A^ = 10. We verified that along the diagonal pi(Q, Q' = Q) = p{Q) 
reproduced the single-particle density |jT6|. The off-diagonal elements of the reduced density 
matrix relate to ODLRO, and it is clear that as A^ increases the off-diagonal elements are 
decreasing in contrast to the diagonal. This is a first indication that ODLRO vanishes for a 
system of hard core bosons in a ID HO in the thermodynamic limit. 



C. ODLRO, natural orbitals and their occupation 

In a macroscopic system, the presence or absence of BEG is determined by the behavior 
of pi{x, x') as |x — x'l — >■ oo. Off-diagonal long-range order is present if the largest eigenvalue 
of pi is macroscopic (proportional to A^), in which case the system exhibits BEG and the 
corresponding eigenfunction, the condensate orbital, plays the role of an order parameter 
|p0| , |2l[| . Although this criterion is not strictly applicable to mesoscopic systems, if the 
largest eigenvalue of pi is much larger than one then it is reasonable to expect that the 
system will exhibit some BEG-like coherence effects. Thus we examine here the spectrum 
of eigenvalues \j and associated eigenf unctions 0j(x) ("natural orbitals") of pi. Although 
natural orbitals are a much-used tool in theoretical chemistry, they have only recently been 
applied to mesoscopic atomic condensates |2^. The relevant eigensystem equation is 



pi{x,x')(j)j{x')dx' = Xj(f)j{x) (14) 

Xj represents the occupation of the orbital (pj, and one has J2j = Numerical evalua- 
tion of the integral (13) by discretization yields a readily-solved matrix eigensystem equation 
which yields accurate numerical results for the largest eigenvalues and associated eigenvec- 
tors. In Fig. ^(a) we show a log-log plot of the fractional occupation of the lowest orbital 
/o = Xq/N versus the total particle number A^ (solid line), along with a best fit power-law 
/o ~ A^~0-^i (dashed line). This is to be contrasted with the case of a spatially uniform 
system of hard core bosons for which /o ~ A^""-^ [|12| . In both cases the fractional occupa- 



tion decreases with increasing A^, and thus do not correspond to true condensates for which 
/o = 1. Nevertheless, the occupation of the lowest orbital may still be large Xq ~ N^-^^^ 
and is larger than the spatially uniform case Aq ~ N^'^, so macroscopic quantum coherence 
effects reminiscent of BEG can still result [p|,p!2|- p!^JT6|JT7[] . 



Figure |(b) shows the distribution of occupations Xj versus orbital number j (the orbitals 
are ordered according to eigenvalue magnitude, the largest eigenvalue being j = 0) for A^ = 2 
(circles), A^ = 6 (stars), and A^ = 10 (squares). This figure shows that as the lowest orbital 
occupation Aq increases with increasing A^ so does the range of significantly occupied higher- 
order orbitals with j > 0. This means that the dominance of the lowest orbital decreases 
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with increasing A^, so singling out (f)o{x) as a macroscopic wave function for the whole system 
becomes more problematic with increasing ||T7|JT6[| . 

The numerically computed lowest orbitals 0o(Q) a^re shown in Fig ^ for a) = 2, b) 
= 6, and c) N = 10, and they show the expected broadening due to many-body repulsion 
as A^ increases. We remark that these lowest orbitals are not simply the square root of the 
corresponding single-particle densities p{Q) as would be the case for a true condensate. 
Figure ^ shows the higher order orbitals (f)j{Q) versus Q for j = 1, 2, 3 and A^ = 10. Although 
the higher-order orbitals differ in detail from the HO orbitals they share the features that 
the orbitals can be chosen real by removal of an overall phase, and that the j^^ orbital has 
j zeros. 



IV. MOMENTUM DISTRIBUTION 

For a spatially uniform system (no trap) the natural orbitals are plane waves, so the 
occupation distribution of the natural orbitals is the same as the momentum distribution. 
Although this is not the case here due to the effect of the harmonic trap potential, the 
momentum distribution is still physically important, so we study it here. In terms of the 
boson annihilation and creation operators in position representation (quantized Bose field 
operators) the one-particle reduced density matrix is 

p,{x,x') = {^B0\i^\x')tij{x)\^B0) (15) 

The momentum distribution function n{k), normalized to J^^n{k)dk = N , is n{k) = 
(\l/so|S^(A;)a(fc)|\I'Bo) where a{k) is the annihilation operator for a boson with momentum 
hk. Then 

n(fc) = (27r)-W dx dx'pi{x,x')e-'^^''-''^ (16) 



oo 



The spectral representation of the density matrix then leads to n{k) = J2j ■^j\pj{k)\'^ where 
the fij are Fourier transforms of the natural orbitals: Pj{k) = (27r)~^/^ 0„(x)e~*'^^(ia; . 
Figure ^ shows the numerically calculated dimensionless momentum spectrum kosc ■ n{K) 
versus normalized momentum k = k/kosc, with kosc = '^T^jxosc^ for a) A^ = 2, b) A^ = 6, 
and c) N = 10. The key features are that the momentum spectrum maintains the sharp 



peaked structure reminiscent of the spatially uniform case [§|,|T2[ for the ID HO, and that 
the peak becomes sharper with increasing atom number A^. This is to be expected since as 
the number of atoms increase the many-body repulsion causes the system to become more 
spatially uniform within the trap interior. 



V. SUMMARY AND CONCLUSIONS 



In summary, we have investigated the ground state properties of a system of hard core 
bosons in a ID HO using the exact many-body wave function obtained using the Fermi-Bose 
mapping theorem. Specifically, we have numerically evaluated the reduced single-particle 
density matrix for the system using Monte-Carlo integration for particle numbers up to 
A^ = 10, and extracted several quantities of physical significance, including the natural 
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orbitals and momentum spectrum. Our main finding is that the lowest orbital occupation 
scales as Aq ~ _/v°-^^, so that the system does not exhibit true BEG, counter to the case of an 



ideal gas in a ID HO [|T9[. Furthermore, this makes the introduction of an order-parameter 
or macroscopic wave function for the whole system more problematic for large A^. We have 
started to seek analytic approaches to derive the observed scaling of the lowest orbital with 
particle with no success so far. We hope that these numerical results may motivate others 
to approach this challenging problem. 

This work was supported by Office of Naval Research grant NOOO 14-99- 1-0806. 
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FIGURES 



FIG. 1. Gray-scale plots of the dimensionless pair distribution function x^g^ ■ D{Qi,Q2) as a 
function of the dimensionless coordinates Qi and Q2, for a) N = 2, h) N = 6, and c) N = 10. 

FIG. 2. Gray-scale plots of the dimensionless reduced density matrix Xosc • PiiQ, Q') as a func- 
tion of the dimensionless coordinates Q and Q', for a) N = 2, h) N = 6, and c) N = 10. 

FIG. 3. Occupation of the natural orbitals: a) fraction of atoms in the lowest orbital /o = Ao/iV 
versus N, and b) Xj versus orbital number j ioi N = 2 (circles), A'' = 6 (stars), and A'' = 10 
(squares). 

FIG. 4. Lowest natural orbitals <Pq{Q) versus normalized coordinate Q for a)N = 2, b) = 6, 
and c) N = 10. 

FIG. 5. Higher-order natural orbits <^j{Q) versus normalized coordinate Q fov N = 10 and a) 
j = 1, b) j = 2, and c) j = 3. 

FIG. 6. Dimensionless momentum distribution kosc " n^K) versus normalized momentum 
K = k/kosc for a)N = 2, b) iV = 6, and c) N = 10. 
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Figure 3 
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Figure 4 
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Figure 5 
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